Assignment 3: Exploring Evictions and Code Violations in Philadelphia¶
Due date: start of class Wednesday, 2/27
In this assignment, we'll explore spatial trends evictions in Philadelphia using data from the Eviction Lab and building code violations using data from OpenDataPhilly.
We'll be exploring the idea that evictions can occur as retaliation against renters for reporting code violations. Spatial correlations between evictions and code violations from the City's Licenses and Inspections department can offer some insight into this question.
A couple of interesting background readings:
1. Explore Eviction Lab Data¶
The Eviction Lab built the first national database for evictions. If you aren't familiar with the project, you can explore their website: https://evictionlab.org/
1.1 Read data using geopandas¶
The first step is to read the eviction data by census tract using geopandas. The data for all of Pennsylvania by census tract can be downloaded in a GeoJSON format using the following url:
https://eviction-lab-data-downloads.s3.amazonaws.com/PA/tracts.geojson
A browser-friendly version of the data is available here: https://data-downloads.evictionlab.org/
# import dependencies
import pandas as pd
import geopandas as gpd
import hvplot.pandas
import cartopy.crs as ccrs
# read the downloaded geojson file
tracts = gpd.read_file("tracts.geojson")
tracts.head()
1.2 Explore and trim the data¶
We will need to trim data to Philadelphia only. Take a look at the data dictionary for the descriptions of the various columns: https://eviction-lab-data-downloads.s3.amazonaws.com/DATA_DICTIONARY.txt
Note: the column names are shortened — see the end of the above file for the abbreviations. The numbers at the end of the columns indicate the years. For example, e-16 is the number of evictions in 2016.
Take a look at the individual columns and trim to census tracts in Philadelphia. (Hint: Philadelphia is both a city and a county).
philly = tracts.loc[tracts['pl'] == 'Philadelphia County, Pennsylvania']
philly.head()
1.3 Transform from wide to tidy format¶
For this assignment, we are interested in the number of evictions by census tract for various years. Right now, each year has it's own column, so it will be easiest to transform to a tidy format.
Use the pd.melt() function to transform the eviction data into tidy format, using the number of evictions from 2003 to 2016.
The tidy data frame should have four columns: GEOID, geometry, a column holding the number of evictions, and a column telling you what the name of the original column was for that value.
Hints:
- You'll want to specify the
GEOIDandgeometrycolumns as theid_vars. This will keep track of the census tract information. - You should specify the names of the columns holding the number of evictions as the
value_vars. - You can generate a list of this column names using [Python's string formatting]:(https://docs.python.org/3.7/library/string.html#format-examples)
value_vars = ['e-{:02d}'.format(x) for x in range(3, 17)]
e_years = ['e-{:02d}'.format(x) for x in range(3, 17)]
print(e_years)
evictions = pd.melt(philly, id_vars=['GEOID', 'geometry'], value_vars=e_years, var_name='year', value_name='evictions_count')
evictions.head()
# change the format of year
evictions['year'] = evictions['year'].apply(lambda x: '20{}'.format(x[-2:]))
evictions.head()
1.4 Plot the total number of evictions per year from 2003 to 2016¶
Use hvplot to plot the total number of evictions from 2003 to 2016. You will first need to perform a group by operation and sum up the total number of evictions for all census tracts, and then use hvplot() to make your plot.
You can use any type of hvplot chart you'd like to show the trend in number of evictions over time.
total = evictions.groupby(['year'])['evictions_count'].sum()
total.head()
total.hvplot(kind='line')
1.5 The number of evictions across Philadelphia¶
Our tidy data frame is still a GeoDataFrame with a geometry column, so we can visualize the number of evictions for all census tracts.
Use hvplot() to generate a choropleth showing the number of evictions for a specified year, with a widget dropdown to select a given year (or variable name, e.g., e-16, e-15, etc).
Hints
- You'll need to use the
groupbykeyword to tellhvplotto make a series of maps, with a widget to select between them. - You will need to specify
dynamic=Falseas a keyword argument to thehvplot()function. - Be sure to specify a
widthandheightthat makes your output map (roughly) square to limit distortions
evictions.hvplot(c='evictions_count', groupby='year', crs=3857, width=500, height=430, dynamic=False)
2. Code Violations in Philadelphia¶
Next, we'll explore data for code violations from the Licenses and Inspections Department of Philadelphia to look for potential correlations with the number of evictions.
We'll be pulling data directly from the CARTO database for the L&I Violations dataset. API information and metadata is available here:
2.1 Count the number of rows in the database¶
First, use the carto package to query the database and count the total number of rows.
from pyrestcli.auth import NoAuthClient
from carto.sql import SQLClient
API_endpoint = "https://phl.carto.com"
sql_client = SQLClient(NoAuthClient(API_endpoint))
query = "SELECT COUNT(*) FROM li_violations"
counts = sql_client.send(query)
counts
2.2 Get a single row of data for column inspection¶
Query the database API and limit your results to a single row. Create a GeoDataFrame from the query results and inspect the columns to identify the column that gives you the date of each code violation.
# select the first item
query = "SELECT * FROM li_violations LIMIT 1"
onerow = sql_client.send(query, format='geojson')
onerow
# convert to geodataframe and check the columns
one = gpd.GeoDataFrame.from_features(onerow, crs={'init': 'epsg:4326'})
one.head()
2.3 Get the data from 2012 to 2016¶
Using the name of the column identified in the previous section, query the database to get all data for years including 2012 through 2016 (inclusive), for 5 years worth of data. Create a GeoDataFrame from the query results.
Notes
- The query operation might take a minute or two to execute.
- The resulting GeoDataFrame should have between 400,000 and 500,000 rows.
query = "SELECT * FROM li_violations WHERE violationdate >= '2012-01-01T00:00:00Z' AND violationdate <= '2016-12-31T23:59:59Z'"
response = sql_client.send(query, format='geojson')
violations = gpd.GeoDataFrame.from_features(response, crs={'init': 'epsg:4326'})
violations.head()
2.4 Remove any violations with missing locations¶
Check if any of the violations have missing data, and if so, trim these rows from the dataset.
Hints
- Take a look at the
isnull()andnotnull()functions. - They can be applied to the
geometrycolumn to test for missing geometries.
# remove missing geometry
violations = violations.loc[violations['geometry'].notnull()]
violations.head()
2.5 Trim to specific violation types¶
There are many different types of code violations (running the nunique() function on the violationdescription column will extract all of the unique ones). More information on different types of violations can be found on the City's website.
Below, I've selected 15 types of violations that deal with property maintenance and licensing issues. We'll focus on these violations. The goal is to see if these kinds of violations are correlated spatially with the number of evictions in a given area.
Use the list of violations given to trim your data set to only include these types.
violation_types = ['INT-PLMBG MAINT FIXTURES-RES',
'INT S-CEILING REPAIR/MAINT SAN',
'PLUMBING SYSTEMS-GENERAL',
'CO DETECTOR NEEDED',
'INTERIOR SURFACES',
'EXT S-ROOF REPAIR',
'ELEC-RECEPTABLE DEFECTIVE-RES',
'INT S-FLOOR REPAIR',
'DRAINAGE-MAIN DRAIN REPAIR-RES',
'DRAINAGE-DOWNSPOUT REPR/REPLC',
'LIGHT FIXTURE DEFECTIVE-RES',
'LICENSE-RES SFD/2FD',
'ELECTRICAL -HAZARD',
'VACANT PROPERTIES-GENERAL',
'INT-PLMBG FIXTURES-RES']
vio_trim = violations[violations['violationdescription'].isin(violation_types)]
vio_trim.head()
2.6 Make a hex bin map¶
The code violation data is point data. We can get a quick look at the geographic distribution using matplotlib and the hexbin() function. Make a hex bin map of the code violations and overlay the census tract outlines.
Hints:
- The eviction data from part 1 was by census tract, so the census tract geometries are available as part of that GeoDataFrame. You can use it to overlay the census tracts on your hex bin map.
- Make sure you convert your GeoDataFrame to a CRS that's better for visualization than plain old 4326.
from matplotlib import pyplot as plt
# project and convert evictions.crs to 3857
evictions.crs = {'init': 'epsg:4326'}
evictions = evictions.to_crs({'init': 'epsg:3857'})
print(evictions.crs)
# convert violations to 3857
vio_trim = vio_trim.to_crs({'init': 'epsg:3857'})
print(vio_trim.crs)
# matplotlib
# create axes
crs = ccrs.epsg('3857')
ax = plt.axes(projection=crs)
# hexbin
hexplot = ax.hexbin(vio_trim.geometry.x, vio_trim.geometry.y, gridsize=50)
# add census geometry
ax.add_geometries(evictions.geometry, crs=crs, facecolor='none', edgecolor='white', linewidth=0.1, alpha=0.3)
# adjust figure size
ax.figure.set_size_inches((13,13))
# add legend
legend = plt.colorbar(hexplot)
legend.ax.tick_params(labelsize=10)
2.7 Spatially join data sets¶
To do a census tract comparison to our eviction data, we need to find which census tract each of the code violations falls into. Use the geopandas.sjoin() function to do just that.
Hints
- You can re-use your eviction data frame, but you will only need the
geometrycolumn (specifying census tract polygons) and theGEOIDcolumn (specifying the name of each census tract). - Make sure both data frames have the same CRS before joining them together!
# select census tracts geometries from evictions
census = evictions.iloc[:,0:2]
census.head()
# spatial join
joined = gpd.sjoin(vio_trim, census, op='within', how='left')
joined.head()
2.8 Calculate the number of violations by type per census tract¶
Next, we'll want to find the number of violations (for each kind) per census tract. You should group the data frame by violation type and census tract name.
The result of this step should be a data frame with three columns: violationdescription, GEOID, and N, where N is the number of violations of that kind in the specified census tract.
Optional: to make prettier plots
Some census tracts won't have any violations, and they won't be included when we do the above calculation. However, there is a trick to set the values for those census tracts to be zero. After you calculate the sizes of each violation/census tract group, you can run:
N = N.unstack(fill_value=0).stack().reset_index(name='N')
where N gives the total size of each of the groups, specified by violation type and census tract name.
See this StackOverflow post for more details.
This part is optional, but will make the resulting maps a bit prettier.
n_violations = joined.groupby(['GEOID', 'violationdescription']).size().unstack(fill_value=0).stack().reset_index(name='violation_count')
n_violations.head()
2.9 Merge with census tracts geometries¶
We now have the number of violations of different types per census tract specified as a regular DataFrame. You can now merge it with the census tract geometries (from your eviction data GeoDataFrame) to create a GeoDataFrame.
Hints
- Use
pandas.merge()and specify theonkeyword to be the column holding census tract names. - Make sure the result of the merge operation is a GeoDataFrame — you will want the GeoDataFrame holding census tract geometries to be the first argument of the
pandas.merge()function.
# The eviction dataset contains duplicated geometries (from 2003 to 2016), so I remove the duplicated ones by selecting census tracts in 2016
census16 = evictions.loc[evictions['year'] == '2016']
census16 = census16.iloc[:,0:2]
print(len(census16)) # there are 384 census tracts in 2016
# merge tracts in 2016 with violations
vio_merged = census16.merge(n_violations, on='GEOID')
vio_merged.head()
2.10 Interactive choropleths for each violation type¶
Now, we can use hvplot() to create an interactive choropleth for each violation type and add a widget to specify different violation types.
Hints
- You'll need to use the
groupbykeyword to tellhvplotto make a series of maps, with a widget to select different violation types. - You will need to specify
dynamic=Falseas a keyword argument to thehvplot()function. - Be sure to specify a
widthandheightthat makes your output map (roughly) square to limit distortions
vio_merged.hvplot(c='violation_count', groupby='violationdescription', crs=3857, width=500, height=430, dynamic=False)
3. A side-by-side comparison¶
From the interactive maps of evictions and violations, you should notice a lot of spatial overlap.
As a final step, we'll make a side-by-side comparison to better show the spatial correlations. This will involve a few steps:
- Trim the data frame plotted in section 1.5 to only include evictions from 2016.
- Trim the data frame plotted in section 2.10 to only include a single violation type (pick whichever one you want!).
- Use
hvplot()to make two interactive choropleth maps, one for the data from step 1. and one for the data in step 2. - Show these two plots side by side (one row and 2 columns) using the syntax for combining charts.
Note: since we selected a single year and violation type, you won't need to use the groupby= keyword here.
# select evictions in 2016
evic16 = evictions.loc[evictions['year'] == '2016']
evic16.head()
# select violation type INTERIOR SURFACES
vio_type = vio_merged.loc[vio_merged['violationdescription'] == 'INTERIOR SURFACES']
vio_type.head()
# plot side by side
evic_plot = evic16.hvplot(c='evictions_count', crs=3857, width=500, height=430, dynamic=False)
vio_plot = vio_type.hvplot(c='violation_count', crs=3857, width=500, height=430, dynamic=False)
(evic_plot + vio_plot).cols(2)
4. Extra Credit¶
Identify the 20 most common types of violations within the time period of 2012 to 2016 and create a set of interactive choropleths similar to what was done in section 2.10.
Use this set of maps to identify 3 types of violations that don't seem to have much spatial overlap with the number of evictions in the City.
# count the number of violations by types and sort descendingly
vio_most = violations.groupby(['violationdescription']).size().reset_index(name='count')
vio_most = vio_most.sort_values(by=['count'], ascending=False)
vio_most.head()
# get the 20 most common types
vio_top20 = vio_most.iloc[0:21]
vio_top20
# filter the violation dataset by 20 most common types
common_types = vio_top20['violationdescription']
common_types
vio_filtered = violations.loc[violations['violationdescription'].isin(common_types)]
vio_filtered.head()
# convert crs to 3857
vio_filtered = vio_filtered.to_crs({'init': 'epsg:3857'})
print(vio_filtered.crs)
# spatial join the filtered violations and census tracts
joined2 = gpd.sjoin(vio_filtered, census16, op='within', how='left')
joined2.head()
# count the number of each violation in each census tract
vio_count = joined2.groupby(['GEOID', 'violationdescription']).size().unstack(fill_value=0).stack().reset_index(name='violation_counts')
vio_count.head()
# merge with the census tract geometries
vio_merged = census16.merge(vio_count, on='GEOID')
vio_merged.head()
# hvplot
vio_merged.hvplot(c='violation_counts', groupby='violationdescription', crs=3857, width=500, height=430, dynamic=False)
# compare violations with evictions
evictions.hvplot(c='evictions_count', groupby='year', crs=3857, width=500, height=430, dynamic=False)
By comparing each of the 20 violation types with the evictions in 2016, it looks like the 3 types that don't seem to be spatially overlapped with evictions are: EXT A-VACANT LOT CLEAN/MAINTAI, LICENSE-RES GENERAL, VIOL C&I MESSAGE.